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Abstract. The lattice Boltzmann method (LBM) has established itself as a valid numerical 
method in computational fluid dynamics. Recently, multiple-relaxation-time LBM has been pro¬ 
posed to simulate anisotropic advect ion-diffusion processes. The governing differential equations 
of advective-diffusive systems are known to satisfy maximum principles, comparison principles, 
the non-negative constraint, and the decay property. In this paper, it will be shown that current 
single- and multiple-relaxation-time lattice Boltzmann methods fail to preserve these mathematical 
properties for transient diffusion and advection-diffusion equations. It will also be shown that the 
discretization of Dirichlet boundary conditions will affect the performance of lattice Boltzmann 
methods in meeting these mathematical principles. A new way of discretizing the Dirichlet bound¬ 
ary conditions is also proposed. Several benchmark problems have been solved to illustrate the 
performance of lattice Boltzmann methods and the effect of discretization of boundary conditions 
with respect to the aforementioned mathematical properties. 


1. INTRODUCTION AND MOTIVATION 

The lattice Boltzmann method (LBM) has gained remarkable popularity as a versatile numer¬ 
ical method for fluid dynamics simulations [Chen and Doolen, 1998]. LBM has its roots in the 
kinetic theory as opposed to the continuum theory. It needs to be emphasized that LBM solves the 
Boltzmann equation instead of solving the continuum field equations. On the other hand, the finite 
element method (FEM) and the finite volume method (FVM) solve the continuum field equations 
directly. Some of the attractive features of LBM are: It can easily handle irregular domains (e.g., 
unstructured pores and fractures in porous media applications), easy to implement even for com¬ 
plicated flow models, and natural to parallelize [Succi, 2001]. Great advances have been made in 
extending LBM to simulate multi-phase flows [Falcucci et ah, 2006], reactive flows [Rienzo et ah, 

2012] , non-linear chemical reactions [Ayodele et ah, 2011], just to name a few. In this paper, we 
limit our scope to LBM-based formulations for advection-diffusion phenomena. 

In the recent years, several key advancements have been made to extend the LBM to simulate 
transport phenomena. To name a few: [Stiebler et ah, 2008], [Shi and Guo, 2009], [Chai and Zhao, 

2013] , [Yoshida and Nagaoka, 2010] and [Huang and Wu, 2014]. Of these works, Yoshida and 

Key words and phrases, lattice Boltzmann method (LBM); meso-scale modeling; multiple-relaxation-time; 
advection-diffusion equations; comparison principle; maximum principle; non-negative constraint; anisotropy; sta¬ 
tistical mechanics. 
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Nagaoka [Yoshida and Nagaoka, 2010], and Huang and Wu [Huang and Wu, 2014] have proposed 
multiple-relaxation-time lattice Boltzmann methods to handle advection-diffusion equations with 
anisotropic diffusivity tensors. 

The governing equations for a transient advection-diffusion system are parabolic partial differen¬ 
tial equations, which possess several important mathematical properties. These properties include 
the maximum principle and the comparison principle [Protter and Weinberger, 1999; Pao, 1993], 
which have crucial implications in modeling physical phenomena. For example, a key consequence 
of the maximum principle in modeling advection-diffusion systems is the non-negative constraint 
of the attendant chemical species. Several factors such as the physical properties of the medium, 
topology of the domain, and the spatial and temporal discretization determine the performance of 
a numerical solution in preserving the discrete versions of the mentioned mathematical properties. 
A discussion on the influence of these factors in the context of the finite element method can be 
found in [Nakshatrala and Valocchi, 2009]. Violations of these mathematical properties can make a 
numerical solution inappropriate for scientific and engineering applications. It has been shown that 
many popular finite element and finite volume formulations for diffusion-type equations violate the 
maximum principle and the non-negative constraint [Liska and Shashkov, 2008; Nakshatrala and 
Valocchi, 2009; Nakshatrala et ah, 2013b]. Recently, numerical methodologies have been proposed 
under the finite element method to satisfy the non-negative constraint and the maximum principle 
by utilizing the underlying variational structure. Since the lattice Boltzmann method does not 
enjoy such a variational basis, these methodologies developed for the FEM cannot be extended to 
the lattice Boltzmann method. 

To the best of our knowledge, the performance of LBM-based formulations for unsteady diffusion- 
type equations with respect to comparison principles, maximum principles, and the non-negative 
constraint has not been documented in the literature. But, such a study is of paramount impor¬ 
tance, as LBM-based formulations are being employed in predictive numerical simulations. We shall, 
therefore, put some popular LBM-based formulations for diffusion-type equations to test, and par¬ 
ticularly show that these formulations violate all the aforementioned mathematical principles and 
physical constraints. The LBM-based formulations of interest in this paper are the non-thermal 
single-relaxation-time LBM for advection with isotropic diffusion, and the multiple-relaxation-time 
methods for anisotropic diffusion proposed in [Yoshida and Nagaoka, 2010] and [Huang and Wu, 
2014]. The mentioned LBM-based formulations have been chosen merely due to their popularity, 
and this paper does not pretend to be exhaustive. 

The remainder of this paper is organized as follows. Section 2 presents the governing equations 
for a transient advective-diffusive system along with a brief discussion on the associated mathe¬ 
matical properties. Section 3 reviews the non-thermal single- and multiple-relaxation-time lattice 
Boltzmann methods for transient diffusion-type equations. In Section 4, several representative test 
problems are presented with a thorough discussion on the numerical results from LBM-based for¬ 
mulations. Section 5 performs a theoretical analysis on the lattice Boltzmann method to obtain a 
simple criterion on the time-step and lattice cell size to satisfy the non-negative constraint. Finally, 
conclusions are drawn in Section 6. On the notational front, a quantity in the continuous setting 
will be denoted by upright font symbols (e.g., u), and a quantity in the discrete setting will be 
denoted by italic font (e.g., u). 
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2. UNSTEADY ANISOTROPIC DIFFUSION-TYPE EQUATIONS 


Consider a bounded open domain Q. We shall denote the boundary of the domain by T = 
where Yl is the set closure of fi. We assume that the boundary T comprises of two parts T N and T D 
such that r N n r D = 0 and T = T N u T D . We denote the part of the boundary on which Dirichlet 
boundary condition is prescribed by T D . Neumann boundary condition is prescribed on T N . A 
spatial point will be denoted by x. The unit outward normal to the boundary is denoted by n(x). 
The time interval of interest will be denoted by [0, T], and the time is denoted by t. For convenience, 
we introduce the following differential operator: 

C [u(x,t)] := + div [v (x, t) u (x, t) - D (x) grad [u (x, t)]] (2.1) 

where u (x, t) is the concentration field. The diffusivity tensor is denoted by D (x), which is assumed 
to be symmetric, positive definite, and bounded above. The advection velocity is denoted by v(x,t), 
which is assumed to be solenoidal. That is, div[v] = 0. The divergence and gradient operators with 
respect to x are, respectively, denoted by div[*] and grad[-]. The initial boundary value problem 
for a transient advection-diffusion system can be written as follows: 


£[u( x ,t)] = g(x,t) 

(x,t) e Yl x (0,T] 

(2.2a) 

u(x,t) = U P (x,t) 

(x,t) er D X [0,T] 

(2.2b) 

(v (x, t) u (x, t) - D (x) grad [u (x, t)]) • n (x) = q p (x, t) 

( X ,t) er N x [o,r] 

(2.2c) 

u(x,t = 0) = Uq (x) 

X £ Yl 

(2.2d) 


where the source/sink is denoted by g(x,t), u p (x,t) is the prescribed concentration, q p (x,t) is the 
prescribed flux, and the initial concentration is denoted by uo(x). It is easy to check that equation 
(2.2a) is a linear parabolic partial differential equation. The initial boundary value problem given by 
equations (2.2a)-(2.2d) satisfies several important mathematical properties, which will be discussed 
next. 

2.1. Mathematical properties. We introduce the following function space: 

f r)n Tii e)^ n 1 

C? (ft x (0,T]) := (u : ft x [0,T] - R | u, —, —, 6 C (ft x (0, T]) J (2.3) 

where C (ft x (0,T]) is the set of all continuous functions defined on Yl x (0,T]. Similarly, one can 
define C(Yl x [0,T]). 

Property 1 (The maximum principle). Let u(x,t) e (12 x (0,T]) n C (ll x [0,T]) be a solu¬ 
tion of the initial boundary value problem (2.2) with dYL = T D . //g(x,t) > 0 then 

min u(x,t) = min min u(x,t), min uo(x) 

(x,t)6^x[0,T] L( X ’ t ) erX [°’ T ] 

Property 2 (The comparison principle). Let ui (x,t) and U 2 (x,t) belong to C\ (Yt x (0,T]) n 
C {Yt x [0,T]). If C [ui (x,t)] > C [u 2 (x, t)] on fix (0,T] ; and u^ (x,t) > (x,t) onYx [0, T] then 

ui (x, t) > U 2 (x, t) on ft x [0, T]. 

Mathematical proofs to the maximum principle and the comparison principle can found in 
[Evans, 1998]. We now show that the non-negative constraint for the concentration can be obtained 
as a consequence of the maximum principle under certain assumptions on the input data. One could 
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(2.4) 






alternatively obtain the non-negative constraint from the comparison principle, which we would not 
present here. 

Property 3 (The non-negative constraint). //g(x,t) >0 in Q, u p (x,£) > 0 onF, and uo(x) > 0 
in Q then 


u(x, t) > 0 in and Vt 


PROOF. Based on the hypothesis, we have the following two results: 

min u(x, t) = min u p (x,t)>0 
(x,t)erx[0,T] (x,t)erx[0,T] 

min un(x) > 0 

xeQ 


These further imply that 


min 


min u(x,t), min un(x) 
(x,t)erx[0 ,T] xeft 


>0 


From the maximum principle, we can conclude that 


min u(x,t)>0 
(x,t)ef2x[0,T] 


which gives the desired result that u(x,t) > 0 in Q and Vt. 


(2.5) 

(2.6) 

(2.7) 

(2.8) 

(2.9) 

□ 


The following integrals will be used in the remainder of this paper: 

J7i(u; fl; t) := f u(x,t) dfl, u;£l;t):= f max [u(x, t), 0] dfl, 

s/ ill s/ n 

J 2 (u;^;t) := f u 2 (x,t)dfl, u.;S7;t) := f (max[u(x,t),0]) 2 dtt (2.10) 

s/ n s/ 


Property 4 (The decay property). If v(x,t) = 0 ; u p (x,t) = 0 on £/&e entire and g(x,t) = 0 
in t/ien 


—J2(u;fi;t)<0 Vt 
dt 

PROOF. Noting that g(x,t) = 0 and v(x, t) = 0 , equation (2.2a) implies the following: 

dn 


f u—— dfl - f u div[D(x)grad[u]] d ft 
Jn dt Jq 


0 


Using Green’s identity and noting that u(x,t) = u p (x,t) = 0 on T, we obtain the following: 

dn 


[ u——■ dfl + [ grad[u] • D(x)grad[u] d ft ■ 
Jn dt Jq 

Noting that D(x) is a positive definite second-order tensor, we have: 

dn 


0 


r dn 

Jn dt 


dft <0 


This further implies that 


_d 

dt 


r d 

/ u 2 (x,t) dfl = — UT 2 (tl; T2; t) < 0 
Jn dt 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 


(2.14) 


(2.15) 


In order to obtain the above result, we have assumed that is independent of t, which is the case 
in this paper. □ 
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In the subsequent sections we will illustrate the performance of some popular LBM-based for¬ 
mulations with respect to the aforementioned mathematical properties in the discrete setting. We 
will also compare the performance of the lattice Boltzmann method with the finite element method 
in this regard. 


3. THE LATTICE BOLTZMANN METHOD 

The lattice Boltzmann method is a numerical method for solving the Boltzmann equation. A 
numerical simulation based on LBM will provide an approximate distribution of particles in the 
discrete configuration-momentum space. Macroscopic quantities such as concentration can then be 
calculated using the obtained discrete distributions. In this paper, we will consider the Bhatnagar- 
Gross-Krook (BGK) model for the collision term: 

|^ + —p-grad[f] + F~ = t(f eq -f) (3.1) 

at m ap A 

where the continuous distribution function is denoted by f = f(p,x,t), the equilibrium distribution 
is denoted by f eq , A is the relaxation time, the macroscopic momentum is denoted by p, and F 
is the external force. We shall use the symbol r for the non-dimensional relaxation time. We 
will now specialize on the LBM-based formulations for isotropic and anisotropic advection-diffusion 
equations. However, for an in-depth discussion on the lattice Boltzmann method, one can consult 
the references [He and Lou, 1997; Succi, 2001]. We first present the single-relaxation-time lattice 
Boltzmann method, and then followed by multiple-relaxation-time methods. 


3.1. The single-relaxation-time lattice Boltzmann method. We will use a uniform and 
structured discretization of the spatial domain 12. The lattice cell size is denoted by Ax, and At is 
the time-step. We will consider a DnQm lattice model, where n is the number of spatial dimensions 
and m is the number of discretized directions for momentum. The lattice models that are employed 
in this paper are shown in Figure 1. The distribution of particles at a lattice node x, at time t 
and along the velocity direction will be denoted by /)(x,£). The evolution of the discretized 
distributions is governed by the following formula: 

fi (x + ei At, t + At) = % (x, t) [translation step] (3.2) 


where the post-collision distribution is denoted by fi and is defined as follows: 

fi(x,t ) = fi(x,t ) - i (fi(x,t) - f- q (x,t)) + WiAtg (x,t) [collision step] (3.3) 

where g(x,t) is the source/sink at a lattice node. The weight associated with the z-th velocity 
direction is denoted by W{. The equilibrium distribution for the z-th velocity direction will be 
denoted by ff q . A popular choice for the equilibrium distribution of an advective-diffusive system 
takes the following form: 

f° q (x,t) = w i u(x,t)^l + V ( X,q (3.4) 


where advection velocity is denoted by v, and c s is the lattice (sound) velocity. The relaxation time 
is given by the following relation: 


D 1 

A tc 2 s + 2 


(3.5) 
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where D is the (isotropic) diffusion coefficient. The concentration at a lattice node can then be 
calculated as: 


m 


u(x,t ) = YjfiiXit) 
2=1 


(3.6) 


We now describe the discretization of boundary conditions. 

(i) Standard method for Dirichlet boundary conditions : Let {j} be the set of unknown distributions 
at a point xT D , and {i} the set of known distributions at that point. We can then write 


fa(x,t) 


W n 


Zp€ {j} W'p 


U J 


'(x,t) 


■ £ /«(*,*) 

> 


Va 6 {j} 


(3.7) 


It will be shown that this standard way of discretizing Dirichlet boundary conditions can 
contribute to the violation of the non-negative constraint. We, therefore, propose a new way 
of discretizing the Dirichlet boundary conditions. 

(ii) Weighted splitting method for Dirichlet boundary conditions : Let x be the spatial coordinates 
of a lattice node on T D . The discretized distributions at each time-level will be assigned as 
follows: 


f a (x,t) = w a u p (x,t) a = 0,-,ra-1 (3.8) 

where u p is the prescribed concentration. 

(iii) Neumann boundary conditions : Let {j} be the set of unknown distributions at a point x e T N . 
We can then discrete the Neumann boundary conditions as follows: 

fa(x,t + At) = fp(x,t) + — q p (x,t + At) Vae{j} (3.9) 

c s Zkzj ek • n (*) 

where /? is the index of a discrete velocity direction such that e a = -eg, and q p is the prescribed 
flux. Note that e a • n (x) < 0 for a e {j}. 


3.2. The multiple-relaxation-time lattice Boltzmann methods. The single-relaxation¬ 
time lattice Boltzmann methods are limited to isotropic diffusion. Multiple-relaxation-time lattice 
Boltzmann methods offer a more suitable framework for simulating anisotropic diffusion, and are 
more easily presentable in the momentum space [Lallemand and Lou, 2000; Li et ah, 2010]. The 
distributions can be transformed to moments using the following linear operation: 

\6i)(x,t) = M\fi)(x,t) (3.10) 

where M is an m x m orthogonal transformation matrix, and | fi)(x,t) is the column vector of 
size m x 1 of the distributions at lattice node x and at time t. The moment corresponding to fi is 
denoted by qi. The BGK lattice Boltzmann equation can then be written as follows: 

I Qi) (x + eiAt,t + At) = 1 6i)(x,t) - S (|^) (x,t) - |e- q ) (x,t)) + AtM\wi)g (x,t) (3.11) 

where S is the m x m matrix of relaxation times. In the case of the single-relaxation-time lattice 
Boltzmann method, the matrix S will be S = I m /r, where I m is the m x m identity matrix, and 
r is defined in equation (3.5). The matrix S need not remain diagonal in the case of anisotropic 
diffusion. In this paper, we shall use the multiple-relaxation-time methods proposed in [Yoshida and 
Nagaoka, 2010; Huang and Wu, 2014]. The treatment of boundary conditions will remain unchanged 
to what was presented earlier (i.e., equations (3.7)-(3.9)). 
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4. REPRESENTATIVE NUMERICAL RESULTS 


In this section, we employ the single- and multiple-relaxation-time lattice Boltzmann methods 
to solve representative diffusion and advection-diffusion problems. For brevity, we shall refer to the 
multiple-relaxation-time method proposed in [Yoshida and Nagaoka, 2010] as the Y-N method, and 
to the multiple-relaxation-time method proposed in [Huang and Wu, 2014] as the H-W method. In 
all the numerical simulations that employ the H-W method, we have taken c\ - 1, C 2 = -2, aq = 8, 
OL2 - ~8 and sq = s\ = ••• = s m = 1 (cf., equations (15)—(IT) in [Huang and Wu, 2014]). In the rest of 
this section, we shall use the following notation: 

^min(^) = min u(x,t) 

xeQ 

U max (t) = max u(x,t ) 

xeQ 

In all the problems that follow, the distributions fi are initialized as follows 

fi(x,t = 0) =WiU 0 (x) (4.3) 


(4.1) 

(4.2) 


4.1. One-dimensional problems. Consider the computational domain to be D = (0,1). The 
diffusion coefficient is taken as D = 1/3, and the advection is neglected. The governing equations 
take the following form: 


f ^ div [Dgrad [u(x, t)]] = g (x, t) 

(x,t) e Q x (0,T] 

(4.4a) 

u(x, t) = U P (x,t) 

(x,t) er D X [0,T] 

(4.4b) 

n • Dgrad [u] = q p (x, t) 

(x,t)er N x[0,T] 

(4.4c) 

u (x, t = 0) = Uq (x) 

xefi 

(4.4d) 


We employ the D1Q 3 lattice model with weights w\ = 1/2, and W 2 = = 1/4. Several different 

cases are considered below. 

4.1.1. Uniform initial conditions. Consider the following initial and boundary conditions: 


Uq (x, t = 0) = 1 Vxefl 

(4.5a) 

IT 

II 

o 

r-h- 

II 

o 

(4.5b) 

u p (x = 1, t) = 0 

(4.5c) 


where t e [0,10 2 ]. According to the maximum principle, the concentration should remain in [0,1] 
at all times. The following conclusions can be drawn from Figures 2 and 3: 

(a) The single-relaxation-time lattice Boltzmann method for diffusion equations violates the maxi¬ 
mum principle. Although the minimum observed nodal concentration is zero, some of the nodal 
concentrations exceeded unity. 

(b) For a given time-step, the violation of the maximum principle can be reduced by reducing 
the lattice cell size Ax. As we will see later, this need not be the case when the diffusion is 
anisotropic. 

4.1.2. Non-uniform initial conditions. Consider the following initial condition: 


u 0 (x,t = 0) 


1 x e [0.4,0.6] 
0 otherwise 


(4.6) 
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Table 1. One-dimensional problem with uniform initial condition and constant 
source: This table compares the numerical solutions under the standard and weighted 
splitting methods for discretizing Dirichlet boundary conditions for several choices 
of Ax and At. In this numerical experiment, one observes that the weighted split¬ 
ting method produces non-negative solutions, and Its accuracy is comparable to the 
standard method of discretization. 


At 

Ax 

l/r 

standard method 

weighted splitting method 

violation 

Error (T) 

violation 

Error (T) 

10“ 3 

10“ 3 

9.995 x 10“ 4 

Yes 

3.14 x 10“ 4 

No 

3.15 x 10“ 4 

icr 4 

icr 3 

1.000 x 10 1 

Yes 

2.98 x 10“ 4 

No 

3.10 x 10“ 4 

icr 5 

icr 3 

9.520 x 10“ 2 

Yes 

2.90 x 10“ 4 

No 

2.94 x 10“ 4 

10“ 6 

10“ 3 

6.667 x 10 -1 

Yes 

2.90 x 10“ 4 

No 

2.90 x 10“ 4 

icr 7 

icr 3 

1.667 x 10° 

Yes 

2.90 x 10“ 4 

No 

2.90 x 10“ 4 


The boundary conditions are: 


u p (x = 1, t) = 0 (4.7a) 

q p (x = 0, t) = 0 (4.7b) 

Even for this problem, the concentration should remain in [0,1]. Maximum and minimum nodal 
concentrations (^min(^) and u max (t)) are shown in Figures 4 and 5. The following conclusions can 
be made from these figures: 

(a) The maximum and minimum bounds can be violated simultaneously. 

(b) Decreasing the time-step while keeping the lattice cell size constant will not alleviate the viola¬ 
tions of the non-negative constraint and the maximum principle. 

(c) Just like the previous problem, the violations vanish with the refinement of the lattice cell size 
for a fixed time-step. 

4.1.3. Uniform initial condition with constant source. We have taken uo(x) = 0 on the entire 
domain D = (0,1), and r D = T with u p (x e T,t) = 0. The time interval of interest is taken as 

T = 1CT 2 . The source is taken as g(x,t) = 1 for x e D and t e (0, T]. The error will be calculated as 

follows: 


Error(t) = — 


n\ 


N 

-UexactOi ^)) 2 


2=1 


(4.8) 


where N is the total number of lattice nodes, and X{ is the spatial coordinate of the z-th lattice 
node. The exact solution to this problem is denoted by u exact (#,£). The obtained numerical results 
are summarized in Table 1, and one can conclude that the weighted splitting method for discretizing 
the Dirichlet boundary conditions produces non-negative nodal concentrations for one-dimensional 
problems. However, the method comes at an expense of marginal decrease in accuracy. 

4.1.4. On comparison principle. Consider the following boundary and initial conditions: 


uo(x, t = 0) = 0 

xefi 

(4.9a) 

u p (x = 0, t) = ul 

te[0,T] 

(4.9b) 

u p (x = 1, t) = 0 

t e [0,T] 

(4.9c) 
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In this numerical experiment, we have taken T - 0.01, and the source term is taken to be zero (i.e., 
g(x, t) = 0). We shall solve the problem using several cases of Ax, At, and ul- Both the standard 
and weighted splitting methods are employed in separate test runs. 

For this problem, if < u^ then the comparison principle implies that u 1 (x,t) < u 2 (x,t) for 
V(x,t) g x [0,T]. A sample result is presented in Figure 6. Several other numerical experiments 
are performed using different choices of Ax and At. Figure 6 and these numerical experiments 
reveal the following conclusions: 

(a) LBM-based formulations do not, in general, respect the comparison principle. 

(b) For one-dimensional problems, refining the lattice cell size Ax for a given time-step At can 
remove the violations of the comparison principle. 

(c) The weighted splitting method does not guarantee the satisfaction of the comparison principle 
even in ID. 


4.2. Two-dimensional problem with anisotropic diffusion on a non-convex domain. 

We now examine the Y-N multiple-relaxation-time method for anisotropic diffusion tensor. The 
computational domain is shown in Figure 7. We have taken L - 1 and T D = r oute r u Fi nner . On 
the inner boundary the prescribed, concentration is prescribed to be unity (i.e., u p (x,t) = 1 for x e 
Tinner)- The flux is prescribed to be zero on the outer boundary (i.e., q p (x,t) = 0 for x e r out er)- 
The anisotropic diffusion tensor is taken as: 


where 


D(x) = Rj DqR <9 


D 0 = 


10 0 
0 10" 3 


(4.10) 


(4.11) 


The orthogonal rotation matrix is denoted by R#, where 6 denotes the angle of rotation. Herein, 
we have taken 6 = tt/ 4. We employed the D2Q5 lattice model. The discrete velocity directions are 
given by 


r [o,o] z = o 

ej = j c [cos ((i - l)7r), sin ((i - 1)tt)] z = 1,2 (4-12) 

[ c [cos ((2i - 5)tt/ 2) , sin ((2i - 5)tt/ 2)] z = 3,4 

where c = Ax/A t. The respective weights are taken as: 

1/3 z = 0 

1/6 i = 1,2,3,4 

The time interval of interest is taken as T - 10~ 2 . We employed the standard method of enforcing 
Dirichlet boundary conditions (see equation (3.7)). Table 2 provides the discretization parameters 
employed in this paper. Figures 8-10 show that the Y-N method violates the non-negative con¬ 
straint. In fact, the obtained minimum concentration is about -0.4, which is a significant violation 
given the fact that the concentration should be between 0 and 1. Another noticeable feature in 
all the cases considered, the minimum concentration converged to a negative value as the time 
progressed. 



(4.13) 
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Table 2. Two-dimensional problem with anisotropic diffusion tensor on a non-convex 
domain: This table provides the minimum concentrations for various discretization 
parameters (i.e., Ax and At). We have taken Ax 2 = At. 


Case 

Ax 

At 

«mm (T) 

1 

1.25 x icr 2 

1.5625 x 10“ 4 

-0.3781 

2 

1.00 X 10“ 2 

1.0000 x 10“ 4 

-0.4072 

3 

5.00 X 10“ 3 

2.5000 x 10“ 5 

-0.4044 


4.3. Two-dimensional problem with anisotropic and heterogeneous diffusion tensor. 

Consider the spatial domain to be 9 = (0,1) x (0,1). We have taken the following anisotropic and 
heterogeneous diffusivity tensor: 


D ( x >y) =e'h + 


ex 2 +y 2 -(l-e)xy 
-(l-e)xy x 2 + ey 2 


(4.14) 


where e « 1 and e' « 1 are arbitrary constants, and I 2 denotes the 2x2 identity matrix. For this 
numerical experiment, we have taken e = 10 -3 and e r = 1CT 10 . The prescribed concentration on the 
entire boundary is taken to be zero. The initial concentration is taken as: 


Uo(x,y) 


1 (x, y) e [0.4,0.6] x [0.4,0.6] 

0 otherwise 


(4.15) 


The time interval of interest is taken as T = 0.025. The H-W method based on the D2Q9 lattice 
model is employed with the lattice velocity c = Ax/At. The discrete momenta are taken as: 


f [0,0] i = 0 

ej = j c [cos ((i - 1 ) 71 / 2 ), sin ((i - 1)tt/2)] i = 1,2,3,4 (4.16) 

( \/2c [cos ((2 i - 9)7r/4), sin ((2 i - 9)tt/4)] i = 5,6,7,8 

with the following weights: 


4/9 

i = 0 


1/9 

i = 1,2,3,4 

(4.17) 

1/36 

* = 5,6,7,8 



The problem is solved using different choices of Ax and At, which are provided in Table 3. This 
table also provides insight on the performance of the H-W method. The spread of the lattice nodes 
that experience violation of non-negative constraint is shown in Figures 11 and 12. The following 
conclusions can be drawn from Figures 11-14 and Table 3: 

(a) The H-W multiple-relaxation-time method violates the non-negative constraint when the diffu¬ 
sion is anisotropic. 

(b) As discussed earlier, the integral J 2 should decrease monotonically with time for pure diffusion 
equations. However, the H-W method does not respect the decay property. Clipping procedure 
can eliminate the violation of the non-negative constraint but does not eliminate the violation 
of the decay property. However, it has been observed that refining the discretization parameters 
(i.e., Ax and At) can improve the performance of numerical solutions with respect to the decay 
property. 

(c) In the case of the integral J7i, noticeable differences appear between the numerical solution with 
negative values and the numerical solution with the negative values clipped. 


10 











Table 3. Two-dimensional problem with anisotropic and heterogeneous diffusion ten¬ 
sor: In this table, the number of nodes that experience violations of the non-negative 
constraint 7V neg , the minimum and maximum observed concentration for various 
choices of Ax and At are shown 

. Note that in all the cases, despite refining the lattice cell size, the non-negative constraint is 
violated. The number of nodes that have negative values for the concentration is denoted by 
7V neg .. The total number of lattice nodes is denoted by N. 


Case 

Ax 

At 

N neg /N x 100 

^min(T^) 

^max(T^) 

i 

5.00 x 10“ 2 

1.00 X 10“ 3 

90/441 x 100 = 20.41% 

-0.0472 

0.5379 

2 

2.50 x 10“ 2 

2.50 x 10“ 4 

762/1681 x 100 = 45.33% 

-0.0311 

0.5576 

3 

1.25 x 10“ 2 

6.25 x 10“ 5 

2867/6561 x 100 = 43.70% 

-0.0193 

0.5768 

4 

1.00 x 10“ 2 

4.00 x 10“ 5 

4207/10201 x 100 = 41.24% 

-0.0144 

0.6188 

5 

5.00 x 10“ 3 

1.00 x 10“ 5 

11894/40401 x 100 = 29.44% 

-0.0068 

0.6021 


4.4. Fast bimolecular reaction. Consider a simple chemical reaction of the form: 

uaA + tibB -> n c C (4-18) 

where A, B and C are the participating chemical species; and tib and nc are their respective 
stoichiometry coefficients. We are interested in the fate of the product C when the time-scale of the 
chemical reaction is much faster than that of the transport processes (i.e., diffusion and advection). 
A detailed description of this mathematical model can be found in [Nakshatrala et ah, 2013a], and 
will not be repeated here. However, the mentioned paper neglected advection in all their numerical 
examples, and the entire paper is devoted to the finite element method. 

4.4.1. Fast bimolecular reaction in a porous medium. We consider the combined effect of trans¬ 
port and chemical reactions in a porous medium. Problems of this type are frequently encountered 
in precipitation studies and ground-water hydrology [Willingham et ah, 2008]. It needs to be men¬ 
tioned that Willingham et al. [Willingham et ah, 2008] have performed numerical modeling but 
employed the finite volume method for the transport problem and LBM for the flow problem. We 
will investigate whether the non-negative constraint will be violated in the numerical simulations of 
such situations under LBM for both flow and transport problems. In the first numerical experiment, 
we shall assume the diffusivity tensor to be isotropic. The advection velocity will be determined by 
solving the incompressible Navier-Stokes equations in the pore structure (see Figure 15). That is, 
obtain v(x,t) and p(x,t) by solving: 

p{~^ + v ’ £ ra( i [ v ] j = ~g ra d [p] + pdiv j^grad [v] + grad [v] T J (4.19a) 

div [v] = 0 (4.19b) 

where v = v(x,t) is the velocity, p = p(x,t) is the pressure, p is the density of the fluid, and p is 
the coefficient of viscosity. The density is taken as p = 1, and the viscosity of the fluid is taken to 
be p = 1CT 2 . The components of the inlet velocity are y^ nlet = 1 and Vy nlet = 0. Diffusion coefficient 
of all the chemical species is taken as D = 10 -2 . A pictorial description of the problem is given in 
Figure 15. The prescribed concentrations for the chemical species A and B are u^ = 1 and = 1. 
The stoichiometry coefficients are taken as ua = 1, ub = 2, and nc - 1. The dimensions of the 
computational domain are L x = 1/2 and L y = 2. The time interval of interest is taken as T = 0.5. 

li 













Table 4. Fast bimolecular reaction in a porous medium: Discretization parameters 
used in the numerical experiment. 


Case 

Ax 

At 

i 

6.25 x 10“ 3 

9.75 x 10“ 4 

2 

5.00 x 10“ 3 

6.25 x 10“ 4 


Table 4 provides the choices of Ax and At employed in the numerical simulation. We employed 
the single-relaxation-time lattice Boltzmann method using the D2Q9 lattice model for both flow 
and transport problems (i.e., see equation (4.16)). The concentration of the chemical species A , 
B and C are shown in Figures 16-17. No negative values for the concentration are observed for 
the reactants A and £>, which is the result of equations (3.6a) and (3.6b) in [Nakshatrala et ah, 
2013a]. However, the concentration of the product C exhibited negative values, which is shown in 
Figure 18. These negative values are very small in comparison with the maximum concentration in 
the domain. We also observed that refining the discretization parameters, At and Ax, can reduce 
the violation of the non-negative constraint, unlike the previous problems with anisotropic diffusion 
tensor. 

4.4.2. Fast bimolecular reaction in an anisotropic and heterogeneous medium. We will consider 
the domain given in Figure 19. Dimensions of the domain are L x = 2 and L y = 1. The advection 
velocity will be derived from the following stream function: 


V’O, y) 


-y - E a k cos ( 

k= 1 ' 


/p fc 7TX 

7T\ . 

_ _ I Cl 

(Qk^y\ 

V L x 

1 ol 
2) 

ii | 

\ L y 1 


(4.20) 


where the parameters are given by 


( Pi,P2,P3 ) = (4.0,5.0,10.0), (qi,q2,q3) = (1-0,5.0,10.0), 


(oq, OL2-, Oifi) = (0.08,0.02,0.01) 


(4.21) 


The components of the advection velocity are calculated as: 


dip(x,y) , x y) 

v x (x,y) =-——, v y (x,y) = +- 


dy 


dx 


The dispersion tensor is taken as: 


D (x, y) = 10 _5 I + p T ||v|| I + (ft - Pr) 


V 0 V 


(4.22) 


(4.23) 


where 0 denotes the tensor product, I is the 2x2 identity tensor, ||-|| is the 2-norm, /3x = 1CT 4 , 
and /?l = 1. The prescribed concentrations are = 1, and the stoichiometry coefficients are 

ua = 1, ub = 2 and nc - 1. The time interval of interest is T = 0.25. We employed the D2Q9 lattice 
model using the H-W method, and obtained the numerical solution for the fate of the product C. 
Discretization parameters for various cases are given in Table 5. Figures 20-22 clearly show that 
the H- W method violated the non-negative constraint, and the violations did not vanish either with 
time or with refinement of the discretization parameters. 


5. A THEORETICAL ANALYSIS 

In this section, we will provide a simple criterion in terms of the discretization parameters to 
satisfy the non-negative constraint for one-dimensional problems. We will also limit our scope to 
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Table 5. Fast bimolecular reaction in anisotropic and heterogeneous medium: Different 
discretization parameters and violation of the non-negative constraint. 


Case 

Ax 

At 

^min (7~) 

^max (7~) 

N ne jN x 100 

i 

5.00 x 10“ 2 

2.50 x 10“ 4 

-0.0209 

0.2704 

275/861 x 100 = 31.94% 

2 

2.50 x 10“ 2 

6.25 x 10“ 5 

-0.0394 

0.3072 

1283/3321 x 100 = 38.63% 

3 

1.25 x 10“ 2 

1.56 x 10“ 6 

-0.0481 

0.3136 

5703/13041 x 100 = 43.73% 


pure diffusion equations (i.e., v(x,t) = 0 ), and dQ = T D . We will restrict the analysis to the D1Q3 
lattice model. We initialize the discrete distributions fi at all lattice nodes as follows: 

fi(x,t = 0) = w;u 0 (» (5.1) 

Since > 0 and uo(x) > 0, we have /^(x,0) > 0. Furthermore, we assume that the Dirichlet 
boundary conditions will be discretized using the weighted splitting method (see equation (3.8)). 
The weighted splitting method guarantees the non-negativity of distributions fi for a lattice node 
on the boundary provided that the prescribed concentration on T D is non-negative. That is, 

if u p (x e r D , t) >0 then fi(x e r D , t) > 0 (5.2) 

So far, we made sure that all the distributions at the previous time-level are non-negative, and the 
discretization of the boundary conditions will not disrupt the non-negativity of distributions. Since 
the distributions at time t are non-negative, equilibrium distributions /) eq will be non-negative in 
the calculation of the collision step at time t + At. That is, 

fi (x,t)>0^u (x , t) = Y J fi (x, t)>0^ (x , t ) = WiU ( x , t) > 0 (5.3) 

i 

If all of these conditions are satisfied, restricting the value of the relaxation time r in equation (3.3) 
can lead to non-negative concentrations at all lattice nodes and for all time-levels. We require that 


1 - 1/r > 0 


(5.4) 


Using equation (3.5) and the above inequality, one can obtain the following condition that ensures 
non-negativity: 


A t> 


Ax 2 

~6D~ 


(5.5) 


That is, if the discretization parameters, At and Ax, satisfy inequality (5.5) then all the distributions 
fi will be non-negative. Non-negativity of all ff s implies the non-negativity of the concentration 
i/(x,£), which stems directly from equation (3.6). 

Note that this result is only valid for one-dimensional pure diffusion equation (i.e., the advection 
velocity is zero) and for D1Q3 lattice model. Furthermore, the above condition does not guarantee 
the preservation of the comparison principles. Deriving similar conditions for more sophisticated 
lattice models in two and three dimensions and for multiple-relaxation-time methods will require a 
more rigorous analysis. 

Another noteworthy point is that we have put stronger conditions on the values of distributions 
fi in order to meet the non-negative constraint for the concentration. In other words, for nodal 
concentrations to be non-negative we made sure that all the distributions are non-negative (i.e., 
fi(x,t) > 0 Vi). However, this condition can be relaxed by allowing some of the fi to be negative, 
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but with an additional constraint that E^o 1 ^ 0. We do not pursue such an approach here, 

but one can consider them in future developments of lattice Boltzmann methods. 

6. DISCUSSION AND CONCLUDING REMARKS 

The maximum and comparison principles are two important mathematical properties of diffusion- 
type equations. The non-negative constraint is an important physical constraint on the concentra¬ 
tion in transport and reactive-transport equations. There are other properties that the solutions to 
diffusion-type equations satisfy under appropriate conditions on the input data; for example, the 
decay property. A main challenge in designing a predictive numerical formulation is to satisfy these 
mathematical principles and the physical constraints in the discrete setting. In this paper, using 
representative numerical examples, we have systematically documented that the current LBM-based 
formulations do not satisfy the maximum principle, the comparison principle, the non-negative con¬ 
straint, and the decay property. We have also shown that the discretization of boundary conditions 
has an affect on the performance of the lattice Boltzmann method in meeting these properties. 
To this end, we proposed a new way of discretizing Dirichlet boundary conditions - the weighted 
splitting method. We then derived a theoretical bound in terms of the time-step and lattice cell 
size that guarantees non-negative values for the concentration under the weighted splitting method 
for one-dimension problems. 

For anisotropic diffusion problems, we considered two representative multiple-relaxation-time 
lattice Boltzmann methods proposed in [Yoshida and Nagaoka, 2010] and [Huang and Wu, 2014]. We 
have shown that these multiple-relaxation-time methods can give unphysical negative values for the 
concentration. It needs to be emphasized that stability conditions for the lattice Boltzmann method 
(i.e., Courant-Fredrichs-Lewy conditions) have satisfied in all our numerical experiments. This 
implies that meeting stability conditions alone does not guarantee the preservation of the mentioned 
mathematical principles in the discrete setting. The main findings of the paper about LBM-based 
formulations can be summarized as follows: 

(a) One-dimensional problems: For a given time-step, one can eliminate the violation of the non¬ 
negative constraint and the maximum principle by refining the lattice cell size. For a given 
lattice cell size, the violation of the non-negative constraint and the maximum principle cannot 
be eliminated by decreasing the time-step. Both these trends are similar to the finite element 
method, which have been reported in [Nakshatrala et ah, 2013b]. 

(b) Critical time-step: Based on a simple theoretical analysis, we obtained the following bound for the 
time-step and lattice cell size to meet the non-negative constraint under LBM for ID problems: 

A 

At ~~6D 

One can obtain exactly the same bound under the single-field Galerkin finite element method 
based on the backward Euler time-stepping scheme for ID problems [Nakshatrala et ah, 2013b]. 
This is an interesting result given the fact that the underlying basis of the lattice Boltzmann 
method (which solves Boltzmann equation to obtain distributions at lattice nodes) is completely 
different from that of the finite element method (which is based on a weak formulation). 

(c) Isotropic vs. anisotropic diffusion: The violations of the non-negative constraint and the maximum 
principle are smaller in magnitude and smaller in terms of spatial extent when the diffusion is 
isotropic. Also, for a given time-step, one can decrease these violations by refining the lattice 
cell size in the case of isotropic diffusion. On the other hand, neither decreasing the time-step 
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nor refining the lattice cell size will eliminate the violation of the non-negative constraint for 
anisotropic diffusion. 

(d) Convex vs. non-convex domains: The magnitudes of the violation of the non-negative constraint 
are larger for non-convex domains. However, it needs to be emphasized that one may have 
violations even on convex domains. 

(e) Comparison principle: The lattice Boltzmann method violates the comparison principle even in a 
simple setting like ID problems or isotropic diffusion. The violation of the comparison principle 
will be more predominant in the case of anisotropic diffusion. 

(f) Decay property: The LBM-based formulations, in general, violate the decay property. 

(g) The lattice Boltzmann method does not posses a variational structure, similar to the one pro¬ 
posed by the finite element method. Due to this reason, the non-negative formulations proposed 
under the finite element method (e.g., [Nakshatrala et ah, 2013b]) cannot be directly extended 
to the lattice Boltzmann method. 

(h) The only procedure that is available to meet the non-negative constraint under the lattice 
Boltzmann method is the clipping procedure, which basically chops off the negative values. 
But, this procedure fixes neither the violation of the decay property nor the violation of the 
comparison principle. Moreover, this method does not have any physical or mathematical basis, 
and it is rather ad hoc. 

One should be wary of violations of the non-negative constraint, and the maximum and com¬ 
parison principles in the numerical simulations using LBM. In the case of isotropic diffusion, the 
authors suggest investigating the occurrence of the mentioned violations, if any of these violations 
occur, they can be significantly reduced by refining the lattice cell size and the time-step in accor¬ 
dance with the CFL condition. However, in the case of anisotropic diffusion, no clear-cut guideline 
for reducing the violations exists. As demonstrated earlier, refining the discretization parameters 
may not improve the numerical solution. A future research direction could be development of 
LBM-based formulations for transient diffusion-type equations (i.e., diffusion, advection-diffusion 
and advection-diffusion-reaction equations) that respect the maximum and comparison principles, 
and meet the non-negative constraint. This paper can serve as a source of benchmark problems for 
such a research endeavor. 
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D1Q3 D2Q5 


D2Q9 


Figure 1. This figure shows the one- and two-dimensional lattice models that have 
been employed in this paper. 
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Figure 2. One-dimensional problem with uniform initial condition: The maximum 
nodal concentration is plotted against time for various time-steps. In all the cases, 
the lattice cell size is taken as Ax = 0.1. For the prescribed data, the maximum prin¬ 
ciple asserts that the concentration should be between 0 and 1 in the entire domain 
and for all times. This figure shows that the single-relaxation-time method violates 
the maximum principle. Moreover, reducing the time-step keeping the lattice cell size 
fixed may not eliminate the violation of the maximum principle. 
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Figure 3. One-dimensional problem with uniform initial condition: The maximum 
nodal concentration is plotted against time for various lattice cell sizes. The time- 
step is arbitrarily chosen as At = 10 -5 for all cases. It can be observed that, for 
a chosen time-step, reducing the lattice cell size can eliminate the violation of the 
maximum principle for one-dimensional problems. 



time 


Figure 4. One-dimensional problem with non-uniform initial condition: The minimum 
nodal concentration is plotted against time for various time-steps. The lattice cell size 
is taken as Ax = 0.1. For the given input data, the maximum principle asserts that 
the concentration should be between 0 and 1 for all times. The numerical solutions 
under the LBM have violated the non-negative constraint. It can be observed that 
reducing the time-step may not eliminate the violation of the non-negative constraint. 
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time 

(b) 


Figure 5. One-dimensional problem with non-uniform initial condition: Maximum and 
minimum nodal concentrations are plotted against the time for various lattice cell 
sizes. The lattice time-step is taken as At = 1CT 5 . The concentration should be 
between 0 and 1. The numerical solutions have clearly violated the maximum (see 
near t = 0) and minimum bounds. However, it is observed that, for a fixed time-step, 
reducing the lattice cell size Ax can decrease the violation. 
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(b) Violation of the comparison principle near x = 0.5 and x = 0.7. 

Figure 6. One-dimensional problem and the comparison principle: The source is taken 
to be zero (i.e., g(x,t) = 0), and the initial concentration is zero in the entire domain. 
The lattice cell size is Ax = 0.1 and the time-step is At = 10 -5 . The prescribed con¬ 
centration at x = 1 is zero, and three different cases are considered for the prescribed 
concentration at x = 1, as indicated in the figure. It can be observed that the LBM 
based on the D1Q3 lattice model violates the comparison principle. 


x-axis 

(a) Violation of the comparison principle near x = 0.3. 
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Figure 7. Two-dimensional problem with anisotropic diffusion in a non-convex do¬ 
main: This figure provides a pictorial description of the test problem. A concentration 
of u p = 1 is prescribed on the inner boundary. 
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Figure 8. Two-dimensional problem with anisotropic diffusion in a non-convex do¬ 
main: This figure shows the concentration profile at time t = 0.01 under the Y-N 
method for the data given by Case 3 in Table 2. The diffusivity tensor is anisotropic 
but spatially homogeneous. Significant negative values for the concentration are ob¬ 
served in the computational domain. 
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(a) Case 1 
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(b) Case 2 
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Figure 9. Two-dimensional problem witR 2 anisotropic diffusion in a non-convex do¬ 
main: The figure shows the regions where the non-negative constraint is violated 
under the Y-N method. 


















Figure 10. Two-dimensional problem with anisotropic diffusion in a non-convex do¬ 
main: The figure shows the variation of minimum value of the concentration with 
respect to time under the Y-N multiple-relaxation-time lattice Boltzmann method. 
It is evident from this numerical experiment that refining the discretization param¬ 
eters, At and Ax, does not alleviate the violation of the non-negative constraint in 
the case of anisotropic diffusion. 


23 

























0.6 

0.5 

0.4 


0.3 

0.2 


0.1 


0 

0 0.5 1 

x axis 

Figure 11. Two-dimensional problem with anisotropic and heterogeneous diffusion ten¬ 
sor: This figure shows the concentration profile at t- 0.025 under the H-W method. 
The parameters for this numerical experiment are given by Case 5 of Table 3. The 
minimum value of u is negative with 'u m in(7~) = -6.8 x 10 -3 . (See the red regions on 
the online color version of this paper.) 
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(a) Case 1 


(b) Case 2 
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(d) Case 4 



(e) Case 5 


Figure 12. Two-dimensional problem with anisotropic and heterogeneous diffusion ten¬ 
sor: This figure shows the regions that have negative values for the concentration 
under the H-W method. 


25 




















0 . 0 ( 


-+o 

£ 


J -0.10® 


rr— : -1- : -1- 

ITl T V.*-*"*-*--*-*-* 

IpT ; y . 

-1 

0 
f ' /j 
.\ / ; 

• i 

j*-*-*?-*-*- 

*- -#- *-* 

.j 

j mwftr'r 

!■■/ i ; \ ! 1 \ i 

ys i • 

\ • k 

x /\ A 

Vv * 

< 

.i.> . 

•r 


< 

i— -• 

Case 1 




i 

4 

►- 4 

Case 2 
i Case 3 

j 



j 


| Case 4 
. Case 5 



-0.20 L 

0.000 0.005 0.010 0.015 0.020 0.025 

time 

(a) Minimum nodal concentration. 
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(b) Total amount of the dispersing chemical species with the 
negative values. 
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(c) Total amount of the dispersing chemical species without the 
negative values. 


Figure 13. Two-dimensional problem with anisotropic and heterogeneous diffusion ten¬ 
sor: The total amounts of the diffusing species under the single-relaxation-time LBM 
and the clipping procedure are plotted a^inst time. The variation of the minimum 
concentration over time is shown in the bottom figure. The parameters in various 
cases are provided in Table 3. 
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(c) Resulting difference from clipping the negative values. 


Figure 14. Two-dimensional problem with anisotropic and heterogeneous diffusion ten¬ 
sor: This figure shows the variation of various integrals defined in equation (2.10) with 
respect to time under the H-W method. The simulation parameters are provided in 
Table 3. It is evident from the figure that the H- W method violated the decay property. 
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Figure 15. Fast bimolecular reaction in a porous medium: This figure provides a 
pictorial description of the test problem. We have taken L x = 1/2 and L y = 2 in the 
numerical experiment. 
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Figure 16. Fast bimolecular reaction in a porous medium: This figure shows the 
concentration profiles of the reactants A and B under the single-relation-time method 
using the D2Q9 lattice model. The diffusivity tensor is isotropic, and no negative 
values for the concentration are observed for the reactants. 
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Figure 17. Fast bimolecular reaction in a porous medium: This figure shows the 
concentration profiles of the product C . 
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Figure 18. Fast bimolecular reaction in a porous medium: This figure shows the 
minimum observed values of the participating chemical species. Concentrations of 
the chemical species A and B did not violate the non-negative constraint. However, 
small negative concentration of the product C are observed. These violations did 
not diminish with time. But, the magnitude of the violation can be decreased by 
refining the discretization parameters. 
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Figure 19. Fast bimolecular reaction in anisotropic and heterogeneous medium: This 
figure provides a pictorial description of the test problem. The reactants A and B 
undergo transport (i.e., both advection and diffusion) and reacts to give product C , 
which in turn gets transported. We have taken L x - 2 and L y = 1 in the numerical 
experiment. 
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Figure 20. Fast bimolecular reaction in anisotropic and heterogeneous medium: This 
figure shows the concentration profile of the product C at t = T = 0.25 under the 
H-W method. We have taken Ax = 1.25 x 10 -2 and At = 1.56 x 10 -6 (see Case 3 in 
Table 5). 
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Figure 21. Fast bimolecular reaction in anisotropic and heterogeneous medium: This 
figure shows the regions where the concentration of the product C is negative at 
t = 0.25 under the H-W method. The violations appear at large number of lattice 
nodes, and are spread across the computational domain. 
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Figure 22. Fast bimolecular reaction in anisotropic and heterogeneous medium: The 

minimum concentration of the product C is plotted against time for various cases 
whose simulation parameters are provided in Table 4. This figure clearly illustrates 
that refining the discretization parameters does not eliminate the violation of the 
non-negative constraint under the H- W method. 


33 



















































